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' The kinetic theory of gases provides methods for calculating Lyapunov ex- 

ponents and other quantities, such as Kolmogorov-Sinai entropies, that char- 
acterize the chaotic behavior of hard-ball gases. Here we illustrate the use of 
these methods for calculating the Kolmogorov-Sinai entropy, and the largest 
I positive Lyapunov exponent, for dilute hard-ball gases in equilibrium. The 

' calculation of the largest Lyapunov exponent makes interesting connections 

■ with the theory of propagation of hydrodynamic fronts. Calculations are also 

Q"^ . presented for the Lyapunov spectrum of dilute, random Lorentz gases in two 

, and three dimensions, which are considerably simpler than the corresponding 

Q>^ I calculations for hard-ball gases. The article concludes with a brief discussion 

' of some interesting open problems. 
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I. INTRODUCTION 

The purpose of this article is to demonstrate that familiar methods from the kinetic the- 
ory of gases can be extended in order to provide good estimates for the chaotic properties of 
dilute, hard-ball gases Q The kinetic theory of gases has a long history, extending over a pe- 
riod of a century and a half, and is responsible for many central insights into, and results for 
the properties of gases, both in and out of thermodynamic equilibrium Strictly speaking, 
there are two familiar versions of kinetic theory, an informal version and a formal version. 
The informal version is based upon very elementary considerations of the collisions suffered 
by molecules in a gas, and upon elementary probabilistic notions regarding the velocity and 
free path distributions of the molecules. In the hands of Maxwell, Boltzmann, and others, 
the informal version of kinetic theory led to such important predictions as the independence 
of the viscosity of a gas on its density at low densities, and to qualitative results for the 
equilibrium thermodynamic properties, the transport coefficients, and the structure of mi- 
croscopic boundary layers in a dilute gas. The more formal theory is also due to Maxwell 
and Boltzmann, and may be said to have had its beginning with the development of the 



We will use the term hard-ball to denote hard core systems in any number of dimensions, rather 
than using the term hard-disk for two dimensional, hard-core systems, etc. 
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Boltzmann transport equation in 1872 At that time Boltzmann obtained, by heuris- 
tic arguments, an equation for the time dependence of the spatial and velocity distribution 
function for particles in the gas. This equation provided a formal foundation for the informal 
methods of kinetic theory. It leads directly to the Maxwell-Boltzmann velocity distribution 
for the gas in equilibrium. For non-equilibrium systems, the Boltzmann equation leads to 
a version of the Second Law of Thermodynamics (the Boltzmann H-theorem), as well as to 
the Navier-Stokes equations of fluid dynamics, with explicit expressions for the transport 
coefficients in terms of the intermolecular potentials governing the interactions between the 
particles in the gas [0. It is not an exaggeration to state that the kinetic theory of gases was 
one of the great successes of nineteenth century physics. Even now, the Boltzmann equation 
remains one of the main cornerstones of our understanding of nonequilibrium processes in 
fluid as well as solid systems, both classical and quantum mechanical. It continues to be a 
subject of investigation in both the mathematical and physical literature, and its predictions 
often serve as a way of distinguishing different molecular models employed to calculate gas 
properties. However, there is still not a rigorous derivation of the Boltzmann equation that 
demonstrates its validity over the long time scales typically used in applications. Never- 
theless, the Boltzmann equation has been generalized to higher density, and in so far as 
they are available, the predictions of the generalized Boltzmann equation are in accord with 
experiments and with numerical simulations of the properties of moderately dense gases . 

In spite of the many successes of the Boltzmann equation it cannot be a strict consequence 
of mechanics, since it is not time reversal invariant, as are the equations of mechanics, and it 
does not exhibit other mechanical phenomena, such as Poincare recurrences Boltzmann 
realized that the equation has to be understood as representing the typical behavior of a gas, 
as sampled from an ensemble of similarly prepared gases, rather than the exact behavior 
of a particular laboratory gas. He also understood that the fluctuations about the typical 
behavior should be very small, and not important for laboratory experiments. To support 
his arguments, he developed the foundations of statistical mechanics, introducing what we 
now call the micro-canonical ensemble. This ensemble is described by giving all systems 
in it a fixed total energy, and then assuming that the probability of finding a system in 
a certain small region on the constant energy surface is proportional to the dynamically 
invariant measure of the small region, given by the Lebesgue measure of the region divided 
by the magnitude of the gradient of the Hamiltonian function at the point of interest on 
the surface. As we know, this micro-canonical ensemble forms the starting point for all 
statistical calculations of the thermodynamic properties of fiuid and many other systems. 

In his attempt to provide a mechanical argument for the effectiveness of the micro- 
canonical ensemble for calculating the thermodynamic properties of fluids, Boltzmann for- 
mulated the ergodic hypothesis, which, in its modern form, states that the time average of 
any Lebesgue-integrable dynamic quantity of an isolated, many particle system approaches 
the ensemble average of the quantity, taken with respect to the micro-canonical ensemble . 
This hypothesis is the subject of several articles in this Volume, but its value for the foun- 
dations of statistical thermodynamics is often questioned ||^. The questionable points can 
be summarized in a few items: (1) The ergodic hypothesis applies to classical systems while 
nature is fundamentally quantum mechanical. (2) Even granting the approximate validity of 
classical mechanics for many purposes, no laboratory system is truly isolated from the rest 
of the universe. Instead, laboratory systems are constantly perturbed by outside influences, 



2 



and these sources of randomness, together with the simple laws of large numbers for systems 
with many degrees of freedom, may be responsible for the utility of the micro-canonical en- 
semble for the calculation of equilibrium properties. (3) Even granting the ability to isolate 
a laboratory system from the rest of the universe, the time it would take for a system's phase 
space trajectory to sample all of the available phase space on the constant energy surface is 
just too long for ergodic behavior to be a physically reasonable explanation for the effective- 
ness of the micro-canonical ensemble. It seems more likely that the equilibrium behavior of 
a system of many particles depends on a number of these factors, including perturbations 
from the environment, the fact that thermodynamic systems have a large number of degrees 
of freedom, and the fact that the physically relevant quantities are projections of phase 
space quantities onto a subspace of a few dimensions. Thus even though the time scale for 
ergodic behavior on the full phase space may be unreasonably long, the projected behavior 
on relevant subspaces may not take very much time for the establishment of equilibrium 
values of the measurable quantities, such as pressure, temperature, etc. Pending the further 
clarification of these and similar issues, it seems fair to say that the complete understanding 
of the reasons for the validity of the micro-canonical ensemble as the basic understructure 
for statistical thermodynamics has yet to be achieved. 

In addition to resolving the various issues needed to complete our understanding of the 
equilibrium behavior of fluids, we would also like to understand the dynamics of the approach 
to thermodynamic equilibrium on as deep a level as possible. Such an understanding would 
enable us to provide a justification of the Boltzmann equation and its many generalizations, 
as well as of the successful use of hydrodynamic or stochastic equations in nonequilibrium 
situations. The counterpart of Boltzmann's ergodic hypothesis for nonequilibrium phenom- 
ena is the assumption that a dynamical system should be mixing, in the sense of Gibbs. That 
is, given some initial distribution of points on the constant energy surface in phase space, in 
a region of positive measure, a system is mixing if the distribution of the points eventually 
becomes uniform over the energy surface, with respect to the invariant measure |]^. If one 
can prove that an isolated dynamical system with a large number of degrees of freedom is 
mixing, then one can show that the phase space distribution function for the system will 
approach an equilibrium distribution at long times, and consequently, quantities averaged 
with respect to this distribution function will approach their equilibrium values. Needless 
to say, the same concerns listed above suggest that the approach to equilibrium may be a 
very complicated affair with a number of possible factors acting in concert or individually 
as circumstances require. One might also ask why a phase space distribution function is 
needed at all, since a laboratory system corresponds to a point in phase space at any given 
time, and not to a distribution of points in phase space. The usual argument given in sta- 
tistical mechanics texts is that it is easier to describe the average behavior of an ensemble of 
points than to solve the complete set of the equations of motion of a single system, and to 
draw conclusions from such a solution. Of course, computer simulations of fluid systems are 
often attempts to solve the equations of motion for an individual system, but they too are 
influenced by noise in the form of round-off errors, and so do not really describe an isolated 
dynamical system, unless one resorts to certain lattice-type models that can be treated by 
integer arithmetic on the computer. Our hope, often unstated, is that the properties we 
explore using the methods of statistical mechanics are somehow typical of the behavior of 
an individual system in the laboratory, even if we know that this cannot be strictly true. 
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We hope, and occasionally can prove, that the deviations from typical behavior are small. 

In any case, it is clearly important to know what role the dynamics of the fluid system 
might play in the approach to equilibrium. Certainly the equilibrium and nonequilibrium 
properties of the system are sensitive to the underlying molecular structure of the fluid and 
to the interactions between the particles of which it is composed PJTO|]. Our experience 
with the Boltzmann equation also assures us that the role of molecular collisions cannot be 
underestimated, even if we are not entirely certain why the Boltzmann equation works for 
an individual laboratory system. Therefore, when faced with a plausible model of a fluid, 
one would like to know if the dynamics of the model is ergodic, mixing, K, or Bernoulli. It 
is of considerable interest to establish the dynamical properties of a large isolated system 
of particles as the starting point for our investigation of the foundations of nonequilibrium 
statistical mechanics, and then to look at the consequences of: (a) external noise on the 
system, and (b) the restriction of our interest to only a small class of functions of the 
dynamical variables needed for physical applications. 

As a step in this direction we show here that kinetic theory can be used to demonstrate 
(not prove) that isolated hard-ball systems are chaotic dynamical systems We will show, 
in fact, that, at low densities at least, hard-ball systems have a positive Kolmogorov-Sinai 
entropy, which we can estimate, and that the largest positive Lyapunov exponent can also 
be estimated. These estimates, in fact, are in good agreement with the results of numerical 
simulations. What these estimates are unable to tell us is whether or not the systems are 
indeed ergodic, mixing, K, or Bernoulli, since we cannot use these kinetic theory techniques 
to show that the phase space consists of only one invariant region and not a countable 
number of them. In fact there is some evidence that if the intermolecular potential is not 
discontinuous but smoother than a hard sphere potential, then there may be some elliptic 
islands of positive measure in the phase space |]T2||. Strictly speaking, the ergodic hypothesis 
is not valid for such potentials. It remains to be seen whether or not this phenomenon is 
ultimately of some importance for statistical mechanics, especially for systems with large 
numbers of particles, and at temperatures where quantum effects may be neglected. 

As a by-product of the analysis given here we will also be able to calculate the largest 
Lyapunov exponents and Kolmogorov-Sinai Entropy for a dilute Lorentz gas with one mov- 
ing moving particle in an array of randomly placed (but non-overlapping) fixed hard ball 
scatterers |]T3|,|T^. This system is much easier to analyze than a gas where all the parti- 
cles are in motion, and was the first type of hard-ball system whose chaotic dynamics were 
studied in detail, either by rigorous methods or by kinetic theory. 

The plan of this article is as follows: In Section II we will set up the equations of 
motion for hard-ball systems that will enable us to analyze the separation of initially close 
(infinitesimally close, actually) trajectories in phase space. The dynamics of the separation 
of trajectories in phase space is, of course, an essential ingredient in analyzing Lyapunov 
exponents and Kolmogorov-Sinai (KS) entropies. In Section III we will apply these results to 
a calculation of the KS entropy of a hard-ball gas using informal kinetic theory rather than 
a formal Boltzmann equation approach. This informal method gives the leading density 
behavior of the KS entropy, but more formal methods are needed to go further. We will 
outline the more formal method based upon an extension of the Boltzmann equation, but we 
will not go too deeply into its solution, since it rapidly becomes very technical. In Section 
IV we outline a method for estimating the largest Lyapunov exponent for a dilute hard-ball 
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gas using a mean field theory based upon the Boltzmann equation |1T5|JT6[] . In Section V 
we apply these methods, both informal and formal, to the calculation of the KS entropy 
and largest Lyapunov exponents for a dilute Lorentz gas with fixed hard-ball scatterers. We 
conclude in Section VI with remarks and a discussion of outstanding open problems. 



II. THE DYNAMICS OF HARD-BALL SYSTEMS 

In this section we will present a method due to Dellago, Posch and Hoover [|17| for 
describing the dynamical behavior of infinitesimally close trajectories in phase space for 
hard-ball systems. We begin with a consideration of a system of N identical hard balls in 
d dimensions, each of mass m, and diameter a. Their positions and velocities are denoted 
by fi and Vi, respectively, where i = 1,2, ... N labels the particles. For simplicity we can 
imagine that the particles are all placed in some cubical volume V = L"^ and that periodic 
boundary conditions are applied at the faces of the cube. The dynamics consists of periods 
of free motion of the particles separated by instantaneous binary collisions between some 
pair of particles. During free motion, the equations of the system are 

n = Vi, Vi = 0. (1) 

At the instant of collision between particles i and j, say, there is an instantaneous change in 
the velocities. It is convenient to write the dynamics in terms of the center of mass motion 

— * — * 

{Rij,Vij) and relative motion {fij, Vij), 

Rij = {ri + rj)/2, Vij 

The change can be described by the equations 

p/ p T// T/ 

J^ij J^ij^ ^ ij Vij, 

r'ij = rij = acr, - = Ms- ■ Vij, 

where the matrix M5- describes a specular reflection on a plane with normal a, i.e., the unit 
vector in the direction from particle j to i at the instant of the collision, 

= 1 - 2aa. (2) 

Nondotted products of vectors are dyadic products and 1 is the identity matrix. In terms 
of the individual particle velocities, the dynamics is described by 

v'i = Vi- {vij ■ d-)a, 

v'. = Vj + {vij ■ a)a. (3) 

These equations (plus the dynamics at the boundaries) are all that one needs in order to 
determine the trajectory of the system in phase space. However, in order to determine 
Lyapunov exponents and other chaotic properties of the system, we need to examine two 
infinitesimally close trajectories, and obtain the equations that govern their rate of separa- 
tion. 



= {vi + Vj)/2, 

= Vi- Vj. 
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Equations for the rate of separation of two phase space trajectories of hard-ball systems 
have been developed by Sinai using differential geometry [^,18|. This leads to an expression 
for the rate of separation in terms of an operator that is expressed as a continuous frac- 
tion. Here we adopt a somewhat different but equivalent approach that is nicely suited to 
kinetic theory calculations. To obtain the equations we need, we consider two infinitesimally 
close phase space points, T (the reference point) and F + 5r (the adjacent point), given by 
{ri,vi,f2,V2, . . .,rN,VN) and by {ri + 6ri,vi + 6vi,r2 + Sr2,V2 + Sv2, . . . ,rN + SrN,VN + SvN), 
respectively. The 2N infinitesimal deviation vectors 6fi, 6vi describe the displacement in 
phase space between the two trajectories. The velocity deviation vectors are not all inde- 
pendent since we will restrict their values by requiring that the total momentum and total 
kinetic energy of the two trajectories be the same. That is. 



N 



N 



6vi = 0, "^Vi- 6vi = 0. 



(4) 



i=l 



i=l 



where, in the energy equation, we have neglected second order terms in the velocity devia- 
tions. We will not use these equations in any serious way since we will be interested in the 
limit of large A^, in which case they are not very important. In the case of the Lorentz gas, to 
be discussed in a later section, the conservation of momentum equation is not relevant, but 
we will require that both of the two trajectories have the same energy, so that the velocity 
deviation vector is orthogonal to the velocity itself. 

Between collisions on the displaced trajectory, the deviations satisfy 



0. 



(5) 



The treatment of the effect of the collisions on the deviation vectors is more complicated. 
One assumes that the two trajectories are so close together that the same sequence of binary 
collisions takes place on each trajectory over arbitrarily long times. Let us suppose that we 
consider a collision between particle i and j on each trajectory. Now, since the trajectories 
are slightly displaced, the collision will take place at slightly different times on each 
trajectory. This slight (in fact infinitesimal) time displacement must be included in the 
analysis of the dynamical behavior of the deviation vectors at a collision. 

The dynamics of the deviations can be considered separately for the center-of-mass co- 
ordinates and the relative coordinates. The center-of mass coordinates behave as in a free 
flight, so that 

= 6Rii. 5VL 



For the relative coordinates, we consider the relative coordinates of two infinitesimally close 
trajectories, rj 



and r*pV*f. 



reference trajectory 


adjacent trajectory 


(collision at t = 0) 


(collision at t = 5t) 


t < 0: 


t < St: 


Vij{t) constant 


v*j{t) constant 


rij{t) = Vijt + aa 


r-^ijit) = v*^it - 6t) + aa* 


t > 0: 


t > 6t: 


Vij{t) = Ms- ■ Vij 


v*^{t) = M,, . v*^ 


'^i'ji^) — "^s- ■ Vijt + aa 


r^^it) = M,, ■ v*^{t - 6t) + aa* 
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The transformations at a collision in terms of the deviation vectors are found with the aid 
of 



where the superscripts and ~ indicate immediately after and immediately before the 
collision, respectively. We should use Eq. (^) in between collisions, i.e., from the last collision 
up to t = 0~. Then we use the collision rule, to be derived in this section, that links and 

to 6ri and 6vi. Eq. (|]) is used again from t = 6t~^, starting with and ^-u-. It is also 
possible to use Eq. (H) from t = 0"*", as if the values of Sflj and Sv'i^j in the above equation 
were valid at t = O''". The difference is an erroneous additional 5t5v'^ for 5rij, when we apply 
Eq. d^) from t = 0^, but this additional term is quadratic in the deviations and therefore 
negligible. In this way, the collision may also be viewed as instantaneous for the deviation 
vectors. 

We can now write 

5rij = a* a — v*j6t — aa = 5aa — Vij6t, 

where 6a = a* — a, and in the last equality we neglect terms quadratic in the deviations. 
Because \a\ = \a*\ = 1, we have 6a ■ a = 0. Taking the inner product with a of the above 
formula gives 

a ■ 6r 



St = -—^^ 6a 



a ■ Vij 



{a ■ Vij)l - Vijcr 
a{a ■ Vij) 



6fij. 



Substitution of these results into 6rlj = a* a — aa — M5- ■ Vij 6t and 6v[j = M5-. ■ v*j — M5- ■ Vij, 
yields 

6^^^ = ■ 6fij, 

6v'ij = -2Qa{i,j) ■ 6fij + ■ 6vij, 

where Q,a{i,3) is the matrix 

a{a ■ Vij) 

In terms of the individual particle deviations, the collision dynamics reads 

= 6fi - {6rij ■ a)a, 



6rj = 6fj + {6fij 



a a. 



6vl = 6vi - {6vij ■ a)a - Qa{i,j) ■ 6fij, 

6v'j = 6vj + {6vij ■ a)a + Q^(z, j) ■ 6fij. (7) 

Eqs. (^, and (|^ are the dynamical equations that govern the time dependence of 
the deviation vectors, {6ri,6vi}. They have to be solved together with the equations for 
the {ri^Vi} in order to have a completely determined system. That is, in order to follow 
the deviation vectors in time, one needs to know when, where, and with what velocities the 
various collisions take place in the gas. 

In the next section, we will use these equations to provide a first estimate of the 
Kolmogorov-Sinai entropy for a dilute gas of hard balls in d dimensions. 
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III. ESTIMATES OF THE KOLMOGOROV-SINAI ENTROPY FOR A DILUTE 

GAS 



We consider hard-ball system from the previous section when the the gas is dilute, i.e. 
na'^ <^ 1, with n the density N/V, and it is in equilibrium with no external forces acting on 
it. Since the hard-ball system is a conservative, Hamiltonian system, one can easily show 
that all nonzero Lyapunov exponents are paired with a corresponding exponent of identical 
magnitude but of opposite sign P JT0|JI9[1 . Obviously there are an equal number of positive 



and negative Lyapunov exponents, and the sum of all the Lyapunov exponents must be 
equal to zero. For the system we consider here, the Kolmogorov- Sinai (KS) entropy is equal 
to the sum of the positive Lyapunov exponents, by Pesin's theorem [^. We will compute 
the KS entropy per particle in the thermodynamic limit, for a dilute hard-ball gas, using 
methods of kinetic theory [pl] , pl[] . 

To carry out this calculation we will use the fact that when an infinitesimally small 
2 A^ci- dimensional volume in phase phase is projected onto the A^rf-dimensional subspace 
corresponding to the velocity directions, the volume of this projection must grow exponen- 
tially with time t, in the long time asymptotic limit, with an exponent that is the sum of all 
of the positive Lyapunov exponents. That is, when we denote this projected volume element 
by SVv(t), as t — >• oo, 

^^^^'^ -expa ^A.^ (8) 




The same result also holds for another small volume element, denoted by 6Vr{t), that is 
the projection onto the iVd-dimensional subspace corresponding to the position directions^. 
The advantage of using an element in velocity space resides in the fact that the velocity 
deviation vectors do not change during the time intervals between collisions in the gas, but 
only at the instants of collisions. We will make use of this fact shortly. 

Our first step will be to obtain a general formula for the KS entropy of a hard-ball system, 
which, in principle, should describe the complete density dependence of this quantity. Then 
we will apply this result to the low density case, using both informal and formal kinetic 
theory methods. 



A. The KS Entropy as an Ensemble Average 

Our object here is to express the KS entropy for a hard-ball system as an equilibrium 
ensemble average of an appropriate microscopic quantity, which in turn can be evaluated by 



^It is worth pointing out that while both 6Vv and 5Vr grow exponentially in time, their combined 
volume, i.e. the original 2A''d-dimensional volume, stays constant. This seemingly paradoxical 
statement can be understood by realizing that almost all projections of the 2iV(i-dimensional volume 
onto A^d-dimensional subpaces will grow exponentially in time with an exponent given by the sum 
of the largest Nd Lyapunov exponents. 
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standard methods of statistical mechanics. To do this we rewrite Eq. 
of a dynamical quantity as 



as a time average 



A»>0 



t 6V^{0) 



lim - / dr — In 'I , 
\dT ""dVM ' 



(9) 



where the angular brackets denote an average over an appropriate equilibrium ensemble to 
be specified further on. Here we have assumed that the hard-ball system under consideration 
is ergodic so that long time averages may be replaced by ensemble averages. Now we can 
use elementary kinetic theory arguments to give a somewhat more explicit form to the 
ensemble average appearing in Eq. (^. Since the volume element in velocity space does not 
change during the time between any two binary collisions, and since the binary collisions 
in a hard-ball gas are instantaneous, the ensemble average of the time derivative may be 
written as 



Ai>0 




da Q{—Vij ■ d-)\vij ■ a\S{fij — aa) In 



5V„ 



(10) 



Here the step function, 9(a;), is equal to unity for x > 0, and zero otherwise. The prime on 
the velocity space volume denotes its value immediately after the collision between particles i 
and j, while the unprimed quantity is its value immediately before collision. In the derivation 
of Eq. ( pUD we consider the rate at which binary collisions take place in the gas, and then 
calculate the change in velocity space volumes at each collision. Thus, we have integrated 
over all allowed values of the collision vector a for a collision between particles i and j with 
relative velocity and included, by means of a delta function, the condition that the two 
particles must be separated by a distance a at collision. The other factors in Eq. (0) take 
into account the rate at which collisions between two particles take place in the gas. A 
more formal derivation in terms of binary collision operators is easy to construct, but not 
necessary for our purpose here |^ . 

Suppose, for the moment, that the iVd- dimensional velocity deviation vector immediately 



after the collision, {6vi^6v2 



. 5v 



Svn) is related to the velocity deviation 



vector immediately before collision through the matrix equation 



f 6vi \ 



A, 



{ 5vi \ 

\ SVn J 



(11) 
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It then would follow that due to the occurrence of the collision 



detAijl, (12) 



and the sum of the positive Lyapunov exponents would be given by 

^ ^^^(^Lzll ^a^-i 1 da&{-v,2-^)\vu-^\6{fi2-aa)\n\detAu\y (13) 

Here we have assumed that the ensemble average is symmetric in the particle indices, and 
chosen a particular pair of particles, denoted by particle indices 1,2. We now must argue 
for the validity of Eq. ([TT| ) and then calculate the determinant of the matrix A12. 

If we examine Eqs. (|^), we see that we can obtain an equation of the form of Eq. (|11|) if 
we can relate the spatial deviation vectors 5rj and Sfj immediately before an {i,j) collision 
to the velocity deviation vectors immediately before collision. Such a relation must be linear 
since we are keeping terms only to first order in the deviations. We therefore make the 
Ansatz that, in general, the spatial deviation vectors are related to the velocity deviation 
vectors through a set of d x d matrices, called radius of curvature (ROC) matrices (even 
though they have the dimension of time), p^;, via 

6m=j:PkrSvi. (14) 

Then some simple algebra shows that 

det A12 = det [M^(l, 2) - Q^(l, 2) ■ (p^ + P22 - P12 - P21)] • (15) 



We can carry the calculation of the KS entropy one step further now by inserting Eq. (|15| 
into the expression for the KS entropy, Eq. (p!3D, to obtain 



^ Aj = ^a"' ^ f dxidx2dp^^dp^2dp2idp22d(y^{.-vi2 ■ (y)\vi2 ■ cf\5{fi2 



aa] X 



In I det [M^(l, 2) - Q^(l, 2) • (p^^ + P22 - P12 - P21)] \^2{xi, X2, Pn, P12, P21, p22)^ (16) 

where Xi = {fi,Vi), dxi = dfidvi, and we have defined a new pair distribution function, 
^2ixi, X2, Pn, P12, P21, P22) by 

^2(2:1, X2, Pn, P12, P21, P22) = - 1) y dx3...dxN Yl' dpijJ='N{xi,X2, xn, Pn, Pnn)- 

(17) 

Here J-'n is the normalized ensemble distribution function for the positions and velocities of 
all of the particles and for all of the ROC matrices. The prime on the product means that 
integrations are not carried out over the four ROC matrices whose indices are 11, 12, 21 and 
22. 

To proceed further, we will need to determine the pair distribution function, JF2, and 
then evaluate the averages needed for the KS entropy. We will do this using both informal 
and formal kinetic theory arguments. We mention here that so far we have made no low 
density approximations, so that Eq. (|1^) is a general expression for the KS entropy of a 
hard-ball system, given the validity of the Ansatz, Eq. (|14D- 
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B. Informal Evaluation of the KS Entropy at Low Densities 



A very simple and useful approximation to the KS entropy [TTI] for a dilute hard-ball 
gas can be obtained by noticing that the spatial deviation vectors for particles 1 and 2, 
immediately before their mutual collision can be written in the form 

Sri = tiSvi + 6fi{Tifi) for i = 1,2. (18) 

Here o is the time of the previous collision of particle i with some other particle in the 
gas, and ti is the time interval between that previous collision and the (1,2) collision. The 
quantity, 5rj(ri,o) is the spatial deviation vector for particle i immediately after the previous 
collision. 

The following argument is correct for dispersing billiard systems such as the Lorentz 
gas, and gives the leading order in the density for the hard-ball gas, system which is only 
semi-dispersing 0. Comparing the two terms on the right-hand side of Eq. (p!8D, one finds 
that their ratio scales, on the average, as the mean free time which is proportional to the 
inverse first power of the gas density. Therefore it seems reasonable that for low densities, 
the terms in Eq. (|T8|) proportional to the ti should be dominant and that the terms 6ri can 
be neglected. If we do this, we can use the approximations 

Pl2 = P21 = 

Pa = til for i = 1,2. (19) 

If we insert this approximation in Eq. (|TB|), we find that 

det Ai2 ^ det [M^(l, 2) - Q^(l, 2)(ti + t^)] . (20) 

The determinant can be evaluated using the explicit expressions for M and Q in Eqs. (^ 
and d^). For two dimensional systems we find 

I det [M^(l, 2) - Q^(l, 2)(ti + t,)] I = 1 + ^''^^^^^^ (21) 

acos0 

where cos0 = |(T ■ 'yi2|/|'5i2|, with — 7r/2 < < 7r/2. For the three dimensional case we find 
I det [M,(l, 2) - Q,(l, 2){t, + t,)] I = 1 + ^^i!l±M(i + eos^ 0) + ( \M^1±M] ' . 



acos( 



(22) 



We can now insert these expressions into the right-hand side of Eq. ( |T^ ) to obtain explicit 
expressions for the KS entropy per particle of dilute hard-ball gases, provided we can specify 



dispersing billiard takes any infinitesimal, parallel pencil of trajectories in configuration space 
before a collision into a defocusing pencil in all directions. A semi-dispersing billiard leaves rays 
still parallel after collision in at least one plane through the pencil. A simple dispersing billiard is 
the outer surface of a sphere, and a simple semi-dispersing billiard is the outer surface of a circular 
cylinder. 
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the pair distribution function, T2- With the approximations made above in Eq. (0), a 
consistent equihbrium approximation for the pair function at low densities is provided by 
the formQ 

J-2(xi,X2,Pn,Pi2,P2i,P22) = ^ Vo (^^1 )</'o (^^2) ^^(^^1 ) ^^(^?2)e-(^("^^)*^+^(^^^)*^)(5(pi2)5 (Pai ) . (23) 

Here v'o('^i) are the normahzed MaxweUian velocity distributions, v{yi) is the collision fre- 
quency for a particle with velocity Uj, 

v{vi) = na'^~^ j dvs J da\vi3 ■ d-\Q{-Vi3 ■ a)(po{v3), (24) 

and we have used the low density expression for the distribution of free flight times, ti 
between collisions. Eq. (^31) is equivalent to the approximation 



J-2(l,2) ^^i(l)^i(2)5(pi2)5(P2i), (25) 

where 

J'i(z) = rupo{vi)i^{vi)exp[-tiu{vi)]. (26) 

It is now a simple matter to calculate h^s/^- To lowest order in density we can expand 
the logarithm of the determinant about the term with the highest power of the time, which 
is linear in time for d = 2, and quadratic in time for d = 3. By the usual scaling arguments, 
the additional terms in the expansion of the logarithm appear to be at least one order higher 
in the density than the terms we keep. Thus, for two dimensional systems 

" ' " " ' " ' " ' COS0|fi2| X 



hxs/N = — — - ~ dvidv2 dti dt2^Poivi)^poiv2) J 

u{vi)u{v2) exp [-{u{vi)ti + iy{v2)t2)] In 



/2 

'\Vi2\{ti+t2) 



acosi 



(27) 



The integral may now be performed, in part, numerically, by using the equilibrium values 
for the velocity dependent collision frequencies. A similar calculation may be done for the 



three dimensional well. In each case we find that [Ol 

hKs/N = AdUd[-\nhd + Bd + 0{nd)]. (28) 
Here is the average collision frequency per particle at equilibrium, 

1/2 = {2TT^/^na)/{f3my/^ and z/3 = {A7r^/V)/{(3my/^, (29) 

where /3 = (ksT)'^, T is the temperature of the gas, and ks is Boltzmann's constant. In 
addition, is a reduced density given by 



^We note that some delta functions are needed to convert p^j from a 2 x 2 matrix to a scalar 
quantity ti, but this is simple to fix, and we do not provide the details here. 
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n2 = na^, and ris = na^vr. 



(30) 



The quantities A^, are numerical factors given by 

A2 = 1/2, and B2 = 0.209, (31) 
A3 = 1, and B3 = 0.562. (32) 

The values for A^ are in excellent agreement with computer simulations by Dellago and 
Posch, but the values for the B^ are too small by factors of 3 or so. These results are 
illustrated in Figs. ([I|) and (|^) together with the results of numerical simulations by Dellago, 
and Posch. 
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FIG. 1. A plot of the results for the KS entropy per particle for a dilute gas of two-dimensional 
hard balls, as a function of the collision frequency ly in units of y/ksT/ (ma). The solid line is the 
result given in Eq. (31), and the data points are the numerical results of Dellago and Posch. Here 
A'^ is the number of particles used in the simulations. The data points are labeled according to 
the computational method, molecular dynamics (MD) or direct-simulation Monte-Carlo (DSMC). 
Also plotted, as the dashed curves, are fits to the data points, to functions of the form (pSD, with 
Ad and B^i as fitting parameters. 
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FIG. 2. A plot of the results for the KS entropy per particle for a dilute gas of three-dimensional 
hard balls, as a function of the collision frequency u in units of ^ ^ksT/ (ma). The solid line is the 
result given in Eq. (32), the data points are results of Dellago and Posch and the dashed curves 
are fits. The notation is the same as in the previous figure. 

The source of the discrepancy between the values of can be attributed to our neglect of 
the spatial deviation vectors in Eq. (pISl). For semi- dispersing billiards, the spatial deviations 
in certain directions remain comparable to the corresponding component of ti6vi. Therefore, 
they may contribute to the first correction to the density logarithm in Eq. (|28|). In the next 
subsection we will pursue this point a bit further. 



C. Toward the Formal Evaluation of the KS Entropy at Low Densities 

As we have seen above, the KS entropy of an equilibrium hard-ball system may be 
calculated if one knows the pair distribution function JF2. The simple approximation to 
this function, Eq. correctly gives the leading order term at low density, but not the 

next term in a density expansion. At low densities, one might expect that the two colliding 
particles, 1, 2 in the above expressions are uncorrelated just before their collision, so that a 
useful approximation to the pair function would still be of the form 

J^2{Xl, X2, Pii, P12, P21, P22) = J^liXl, pn)^l(a;2, P22)^(P12)^(P21)> (33) 
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but we should use a better approximation to the one particle distribution function, jFi, than 
that used above. More specifically, we consider the set of reduced distribution functions, 
jFi, JF2, which satisfy a set of BBGKY hierarchy equations, and then use this set to find 
a good, low density approximation to the equation for in much the same way as the 
BBGKY hierarchy is used to derive the standard Boltzmann transport equation and its 



higher density corrections [^. However as we now consider a set of extended distribution 
functions which include the p matrices as variables in addition to the positions and velocities 
of the particles, the hierarchy equations must be extended to include these new variables 
as well. The full hierarchy technique requires the use of somewhat cumbersome cluster 
expansion methods, and has been outlined elsewhere |]21|. Nevertheless, the low density 



approximation to the equation for T\ is physically plausible since it is a direct extension of 
the Boltzmann equation to include the additional ROC matrices as variables. This extended 
Boltzmann equation is given by 



J'l(Xi,Pii,t) rfX2rfpi2C?p21^P22^-(l>2)J^i(Xi,Pii,t)j^l(x2,p225i^) >< 

<^(Pl2)5(P2l)- (34) 



Here we have included a possible time dependence of the distribution functions, the operator 
£0(1) is the free streaming operator which accounts for the change in T\ due to the free 
motion of the particle. It is given by 

Here we assume that there are no external forces acting on the particles, so that in the course 
of free motion for a time interval, (r, t + r), the position of the particle changes according 
to ri(t + r) = ri(r) +tvi(r), and the ROC matrix, according to P;^(t + r) = Pi(t) The 
derivatives with respect to the position variable and with respect to the diagonal components 
of refiect this free fiight motion. The right-hand side of the extended Boltzmann equation, 
Eq. (P^D, expresses the change in T\ due to binary collisions, and we have supposed that 
the two colliding particles are uncorrelated before their collision, as is normally done in 
the derivations of the Boltzmann equation. In this term we have defined a binary collision 
operator by 



r_fi,2i 



da\v\2 ■ O"|0(^^i2 ■ 0") \^{fx2 - acr) J dp[^dp[2dp2idp22 

X n 5(p.,-p.,(p'n,...,py)n(i,2) 



-5{fi2 + aa) . (36) 

Here Vg. is a substitution operator that replaces ROC matrices to its right by the corre- 
sponding primed matrices, and velocities to its right by restituting velocities, namely those 
which for a given a produce vi,V2 after a binary collision. The delta functions in the ROC 
matrices require that the primed ROC matrices be restituting ones, namely, those which 
produce the unprimed matrices after the 1, 2 binary collision. 
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The restituting values of the ROC matrices can be found by means of Eqs. (H) and 
and the Ansatz, Eq. (0). To do this we consider only the ROC matrices involving the two 
colliding particles, and ignore the other ROC matrices involving other particles, such as p^a? 
etc. In the equations below we will denote with primes the values of restituting quantities, 
i.e., the quantities before collision will be denoted with primes, and those after collisions 
without primes. Then before collision 

= p'u ■ M + P'u ■ ^^2^ 

6r^ = p'22 ■ + P21 ■ K? (37) 

and it follows that 

6fi, = ■ 6vi, + p', ■ 5V;„ (38) 

where 

Pa= 2 + + ^21 + P22) > 

Pb=\ (Pll - P12 + P2I - P22) ? 
Pc=\ (PU - Pl2 - P2I + P22) ? 

Pd = (Pii + P12 - P21 - P22) • (39) 

Now one writes the collision equations for the spatial deviations of the center of mass and 
relative coordinates in terms of the variables before (primed) and after collision (unprimed). 

Pa ■ ^Vl2 + Pb ■ ^^12 = P'a ' ^^12 + P'b " ^^12, 

p, ■ 6v^2 + P, ■ SV^2 = M, ■ [p^ ■ 6^f,, + p', ■ 6V;,], (40) 
as well as the collision equation for the deviation of the relative and center of mass velocities, 

5vu = M, ■ 51^12 - 2Q', ■ [p^ ■ K2 + P'd ■ 5V;,], 

SVu = SV;^. (41) 

Then we use the fact that the deviations of the relative and center of mass velocities before 
collision are independent of each other, and that there are no orthogonalization or nor- 
malization constraints on the components of these vectors that would make some of them 
dependent on the other^. Therefore, upon inserting Eq. (^Tj) into Eq. (|40|) , and comparing 
coefficients of 6V(2 and of Sv'12, we obtain 



^The constraints mentioned earlier, Eq. (^), on the deviations of the total momentum and energy 
of the A'^-body system do not have any significant effect on the dynamics of two particles if ^ 2. 
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Pa 


-2p, 


q; 


p'd 


= P'a^ 






-2p, 


q:. 


P'c 


= pI 






-2Pe 


q:. 


P'c 


= M^ 


p'c. 


Pd 


-2Pe 


q; 


p'd 


= M^ 


p'd- 



(42) 

In order to complete the specification of tfie restituting ROC matrices, we use tlie fact 
tliat the particles are taken to be dynamically uncorrelated before collision so that p'^2 — 
p'21 = 0, and we can also express the matrix in terms of the unprimed relative velocity 
after collision: = Qa{v[2) = —Qai^i'^)- Thus, the primed variables before collision can 
be expressed in terms of the unprimed variables after collision, provided the Eqs. (^) for 
the p' matrices can be solved. It is important to point out that, as we will see in Section V, 
the solution of these equations is very different for dispersive billiards such as encountered in 
the Lorentz gas, and semi- dispersive billiards as occur in the hard-ball gas where all of the 
particles move. The difference resides in the fact that for dispersive billiards the elements 
of ROC matrices after collision are of the order of a/v and therefore much smaller than the 
typical matrix element before collision, which is on the order of the mean free time. For 
semi-dispersing billiards this is no longer true, and some of the matrix elements remain large 
after collision as well. This is equivalent to the statement that the spatial deviation vectors 
immediately after collision are not negligible for semi-dispersing billiard systems. 

It is an elementary matter now to write expressions for the low density KS entropy in 
terms of the distribution functions and the elements of the ROC matrices. One need only 
replace JF2 in Eq. (|16D by J-'i{xi, Pii)J-'i{x2, P22)'^(Pi2)^(P2i)) use the leading order term 
in the expression for the determinants, namely 

I det [M,(l, 2) - Q,(l, 2)(pn + P22)] I ^ ^J^^l^, (43) 

acoscp 

for two dimensional systems. Here the matrix element p±± is defined by 

P±± = ^{Vl2L ■ (Pll + P22) ■ ^12±)- (44) 

Here the subscript ± on the two dimensional unit vector vul denotes a unit vector in a 
direction perpendicular to Vu- For the three dimensional system, we obtain 

I det [M,(l, 2) - Q,(l, 2)(p,, + P22)] I ^ [pWp'A - P±±pI\]- (45) 

Here 

pI± = ^(^l2± ■ (Pii + P22) ■ ^{21.), (46) 

where the unit vectors t)i2, t'i2±? ^i2± form an orthonormal coordinate basis in three dimen- 
sions. 

It now remains to solve the extended Boltzmann equation for JF^. This appears to be 
somewhat difficult due to the complications inherent in the restituting, or gain term. If 
one ignores this term entirely, for whatever reason, the solution becomes simple. In fact 



one recovers the approximate form for J-'i given by Eq. (26), and therefore one recovers the 
same expressions for the KS entropy as given by Eqs. ([25|-|5^ pT| . To do better, one has 
to include the restituting term in the equation for J-'i. How to (approximately) solve the 
extended Boltzmann equation then, is still under investigation. 
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IV. THE LARGEST LYAPUNOV EXPONENT FOR THE HARD-BALL SYSTEM 

AT LOW DENSITY 



In this section, we will obtain an estimate for the largest Lyapunov exponent Xmax- This 
exponent determines the asymptotic growth of all the deviation vectors: 6vi{t) , 6ri{t) ~ 
gAmaii ^jj] calculate it in the low density regime, -hii <^ 1, when a Boltzmann equation 
gives an appropriate description of the behavior of the one-particle distribution function. 
The results will be checked against simulations. 



A. Kinetic Theory Approach 

At low density, two colliding particles, i and j, have spent a long time in free flight before 
the collision. Therefore, in Eq. (|18D, just before collision we keep only the dominant term 
involving the free flight time ti. We saw in the previous section that for the KS entropy 
this approximation turned out to be unsatisfactory for describing the behavior of the ROC 
matrices. For the leading behavior of the velocity deviations, however, it will be a valid 
approximation. 

We insert this in Eq. (0) and keep only the leading terms in the free flight times, 

Sr^ ^ -6rj ^ -iM^iUSvi - tj6vj), 

5v'^ ^ ^ q,{tM - tj6v,). (47) 



Let us denote the typical velocity of a particle in the gas by Vq = ^JkBT/m. One sees that 
Sv'Jvq and Sfja are of the same order order of magnitude, which justifies the approximation 
that 6fi -C ti6vi. Notice we have eliminated the position deviation vectors, but now we need 
to know the free flight times between collisions. To compare different contributions in 
Eq. (p7D, we want to know the order of magnitude of the velocity deviation vectors. We 
define clock values ki to represent their order of magnitude measured in inverse powers of 
rid-, 

5vi = 5vo{ — ] Ci, (48) 



where Cj is a unit vector and 6vq is a fixed infinitesimal velocity. To obtain the largest 
Lyapunov exponent, it is enough to know the time evolution of these clock values. They 
will grow linearly with time, 

ki ~ wudt- 

The average collision frequency I'd sets a time scale proportional to the density. It is extracted 
so that the clock speed w can be interpreted as the average increase of a clock value in a 
collision. We see from Eq. (^) that the clock speed is related to the largest Lyapunov 
exponent via Xmax = —wud^^nd- 

After the collision, the two particles have equal clock value, given by 



j^, _ j^, _ In \Qa{i,j){tiSvi - tj6vj)/6vo 



(49) 

max(A;i, kj) + 1 + C (1/ Inn^) . (50) 
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To see how Eq. (^0]) follows from Eq. (|49|) consider the following. If one of the clock values, 
say ki, is larger than the other one, the corresponding velocity deviation 6vi ~ n^'^' is very 
much larger than 6vj ~ ■'. We will consider the case that the two clock values differ by 
less than one in a moment. We say particle i is the dominant particle. Only the dominant 
term inside the logarithm in the denominator of Eq. (|49|) has to be taken into account, as 
the other one is much smaller. The free flight times scale, for low densities, inversely with 
the density hd, whereas the terms Qa{hj) and vq do not contain any density dependence. 
To see the density dependence more explicitly, we scale the free flight time as = hdti, so 
that Si is density independent. The denominator in Eq. (^) is then seen to be the sum of a 
term proportional to — In and a density independent term which fluctuates from collision 
to collision: 

lia\Qa{i,j)ti6vi/6vo\ = \n\Qa{i, j)si{l/nd)'''^^ei\ = {ki + In ha) + \n\Qa {i, j)siei\, 



which leads to Eq. (|50D . The fluctuating corrections to the addition of 1 of the dominant clock 
value scale with density as 1/ Inn^. Finally, let us consider what happens when \ki — kj\ < 1. 
In that case, it is not necessarily true that one of the two terms in the denominator in 
Eq. dominates. To show that Eq. (^Up then still holds, i.e. /c- = kj = max{ki, kj) + 1 
plus corrections which scale with the density as l/lnn^, we consider the extreme case that 
ki = kj exactly: 

^ \ki+li 



In \Qa{i J) {USvi -tjSvj/5vo\ = \n\Qa{i,j){l/nd) {siCi - SjCj^ 



{ki + l)(-lnnrf) + In \Qa{iJ){siei - SjC 



where Sj = hdtj. Again, the second fluctuating term is density independent, and Eq. (|50D 
follows. 

The clock values have a well-defined dynamics in the limit fid 0, which gives the 
leading behavior of w. For this leading behavior, we can neglect the 0{1/ Inhd) terms in 
Eq. (li: 

fc- = k'j = m.ax{ki, kj) + 1. (51) 

In some collisions the neglected terms may be large, but the lower the density, the rarer 
such collisions become. If these terms were kept, they would give rise to terms of the order 
of \n~^{nd) and higher in w, so we would get an expansion of the form: 

w = Wo + Wi \n~^{nd) + W2 \n~'^{nd) + 

The coefficients Wq, Wi, etc., in fact are still density dependent due to the neglected terms 
in Eq. {^7\j, and to correlations between collisions. The former will give rise to powers of 
the density, the latter may also give rise to a nonanalytic dependence on the density ^3 



Such terms however are at least accompanied by one factor of the density. That means 
that for "hd — >■ 0, the limiting values of the coefficients give the asymptotic behavior of w. 
Corresponding to the expansion of w, the Lyapunov exponent has an expansion of the form 

Amax = i^d[-wo\nhd -wi- W2ln~'^{hd) + ...]. (52) 

We will calculate the first term in the density expansion of Xmax (which gives the leading 
behavior) by calculating the leading clock speed Wq. As we already argued, it is enough to 
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use the dynamics in Eq. (|5T|), for which we can restrict ourselves to integer clock values. 
The calculation will be based on an extended Boltzmann equation for the one particle 
distribution function f{k, v, t) of clock values and velocities. Colliding particles are assumed 
to be uncorrelated before the collision, i.e., the probability that a particle with clock value 
ki and velocity Vi and a particle with clock value kj and velocity Vj collide is proportional to 
f{ki,Vi,t) f{kj,Vj,t); this is the Stosszahlansatz. For that collision we still need to specify 
a and this vector will, in a collision with given Vij, be drawn from a probability distribution 
proportional to \a ■ Vij\. We demand that a ■ Vij be negative, so that the particles are 
approaching each other. The gas will be uniform in equilibrium, therefore we neglect any 
spatial dependence. Considering gains and losses, we can construct the following extended 
Boltzmann equation for the time evolution of /: 



dt 

+f{k~l,v'2) E /(/,t;i) + /(fc-l,^i)/(fc-l,^2 

Z=— oo 



l=—oo 

(53) 



where in the integral, primed and unprimed quantities denote values before, respectively 
after collision. i^{vi) is the velocity dependent collision frequency, given by Eq. (^4|), where 
we take the velocities to have their equilibrium distribution ip^. In principle, we could 
also allow for an initially nonequilibrium or nonstationary velocity distribution, but this 
will not influence the asymptotic clock speed, as the velocity distribution function will be 
stationary in due course, whereas the clock values keep on growing indefinitely. The extended 
Boltzmann equation can be simplified using cumulatives, defined by 

1 k 

Because of this definition, C{k, v) will tend to 1 as tends to infinity for all v. In terms of 
the cumulatives, the Boltzmann equation reads 

— = -u{vi)C{k,vi) 

+ J dv2j rfae(a-{;i2)na''-^|a-t;i2|<^o(^2)C(A;-l,tri)C(A;-l,tr2). (54) 

This is the appropriate equation to calculate Wq from. The assumptions made in its deriva- 
tion were that colliding particles are uncorrelated, that the clock values increase according 
to (pT|), that spatial fluctuations do not play any significant role and that the velocity dis- 
tribution has reached its equilibrium form. The equation is similar to that found in earlier 
work ||T5| , [T6[| , where we had not accounted for the velocity dependence. We will follow the 
same approach used there, to find a better estimate for wq, and thus for the leading behavior 
of the largest Lyapunov exponent for low density. 

The Boltzmann equation is of the same type as equations encountered in front propaga- 



tion |2^. In this analogy, C = 1 is an unstable phase, that lives at the higher clock values. 



and C = is a stable phase, that propagates from the lower clock values to the higher 
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ones. Because the 6vi grow exponentially, this propagation means that the clock values are 
increasing linearly in time: 



C{k,v;t) = F{k- WQUdt, v) . (55) 

We assume that the equation is of the pulled front type, meaning that the equation linearized 
around the leading edge {F = 1) sets a critical velocity of the front. Fronts exist for any 
velocity above this critical one. If the initial conditions of C are sufficiently steep, then 
a front will develop with this critical velocity. In our case, as an initial condition, almost 
any (single) ST will do to see an exponential growth with the largest Lyapunov exponent, 
as almost any ST will have some component along the fastest expanding direction in phase 
space. Because of the finite number of particles, the initial distribution of clock values 
corresponding to a single ST, will have a finite support. This is as steep an initial condition 
as one can get, so the critical velocity is the clock speed Wq that determines the Lyapunov 
exponent. 

We will now show how the linearized equation sets a critical velocity, and then we will 
determine it for the hard-ball gas in two dimensions. We insert the Ansatz (|55|) into the 
Boltzmann equation and concentrate on the leading edge. That is, we write F{k,v) = 
1 — A{k,v) and get, to linear order in A, 

-woUd = -u{vi)A{k,vi) (56) 

+ Jdv2j da e {a ■ vi2)na'^'^\a ■ vi2\M^2) [A{k ~ l,v^) + A{k - l^tf^)] . 
We now take a linear superposition of exponentially decreasing functions in k: 

Aik,v) = Y.Aiv)e-^'\ (57) 



For this to be a solution of Eq. (|5^), only certain characteristic values 7^ and corresponding 
characteristic functions Ai{v) are allowed, determined by the characteristic equation: 

Wo-fi'dA{vi) = -u{vi)A{vi) + I dv2 / da<c){a ■ Vi2)na'^'''\a ■ Vi2\^Po{v2) [A{v'^) + ^({^2)] • 



which is found by inserting A('y) exp(— 7A;) into the linearized equation, Eq. (|56D . The 
characteristic equation is non-linear in 7. To deal with it, we rewrite this equation first as 
a generalized eigenvalue problem: 

A W^,^ A = LA, (58) 
where A = e~'^ is the eigenvalue and the operators are defined by 
{W^,^A) (iTi) = {ytv{vi) + wo-i)A{vi), 

{LA){vi) = z/ji / dv2 [ dae{a-vi2)na'^-'\a-vi2\Mv2) [Aitr,) + A{v'^)] . 



For solving Eq. (0), first the whole spectrum {A^lwy)} of the generalized eigenvalue prob- 
lem (|58|) should be determined, with the corresponding eigenf unctions {A^^v] WQ'y)} . The 
eigenvalues should also equal e^. So next, for every A^, the solutions for 7 of 
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for fixed wq, have to be found, and they are characteristic values. For each term in Eq. 
the 7j is one of these solutions 7 with characteristic function Ai{v) = y4fc(w;wo7j). The 
superposition (|57D , with arbitrary coefficients, then is the general solution of the linearized 
equation (|56|), for a given clock speed Wq. 

However, we don't need the full solution of Eq. (^). For the leading edge, k 00, 
only the term with the slowest decay in Eq. (|57D survives. The characteristic value with 
the smallest real part, we call this 70, corresponds to the slowest decay. As the distribution 
function / should be positive, F should be an increasing function and A should be mono- 
tonically decreasing, meaning that 70 cannot have a imaginary part. The requirement that 
7o be real is essential: it will be the determining factor for the critical velocity wq, as follows. 
The smallest characteristic value 70 is a function of the clock speed wq. The inverse of this 
function, ^0(70), has a minimum. That minimum is the critical value, because, for clock 
speeds Wq below this minimum, there are no real 70. 

The characteristic value with the smallest real part, was required to be real, so we only 
need to consider the eigenvalue problem ( |55| ) for real wq'j. In that case, the operators are 
self-adjoint with respect to the inner product 

{A\B) = J dvipo{v)A{v)B{v). (59) 

The self-adjointness means that all the eigenvalues are real. As we want the smallest real 
7, we are interested in the largest eigenvalue, which we will call Aq. The self-adjointness of 
the operators, together with the fact that W^o-r is a positive operator, can be used to derive 
a maximum principle for this eigenvalue, 

A„K7)=max ^j^^^iy^^ . (60) 

The corresponding eigenvector Aq{v) is the A{v) for which the expression takes on its max- 
imum value. Once we have Ao(wo7); 7o is the smallest real solution of 

A0K70) =e-^°. (61) 



Notice that 70 then indeed still depends on the value of Wq in the Ansatz (55) 



As in previous work [[15| , p!6| , the critical velocity is determined as the smallest Wq possible 



when 7o is real and positive, i.e., the minimum of Wo(7o)- This minimum is obtained from 
Eq. ( |6TD taking the derivative with respect to 70 and using ti'o(7o) = 0. This gives 

woAo(i(7o7o) = -6"^°. 

This condition together with Eq. (pi]) can be captured concisely by saying that the critical 
velocity is the minimum of the following function w over real positive values of Wq'Jq: 

= -lnA„(u>„o.) 

Once the critical clock speed Wq has been obtained, the characteristic function Aq{v), 
determined by Eq. (0), gives the velocity distribution in the head of the front. More 
precisely. 
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(63) 



is the velocity distribution of the particles with clock values k larger than some k^, for k^ 
tending to infinity. 

We have no exact expression for the eigenvalue Aq, but the variational principle will allow 
us to calculate approximations to it, as follows. Aq was given in Eq. (^D|) as the maximum 
over all functions A{v). So inserting any function A{v) would give a lower bound on Aq. 
Using this lower bound in Eq. (^) will then give a lower bound on wq. We take a fixed form 
of A{v; bo, bi, . . .), depending on variational parameters and determine the maximum 
in Eq. (^) over these parameters. This gives an approximation to Aq and Aq{v), which 
can be systematically improved, and checked for convergence, by including more variational 
parameters. 

We will now apply this procedure explicitly to the two dimensional case. The extraction 
of the factor Ud in the definition of L and W^o'y removes the density from the equations. We 
can get rid of the temperature dependence by rescaling the velocities by the typical velocity 
Vq-. V v/vq. The operator W^^-y then becomes 



W^,yA){v^) 



Wo7 + 



2v^ 



27T 



-2\v2i\ A{vi) 



and L is 



{LA){v,] 



2^ 



dvo 



27T 



da\a- V2i\ [A{v'^) + A^v'^)] ■ 



As our variational eigenfunction we choose: A{v) = a + 6o|'^P + ^il'S'l^ + b2Vx + b^Vy. a is 
not a variational parameter, but is set by the normalization requirement, (A |1) = 1. The 
operator L contains no preferred direction, so the eigenfunctions are expected not to contain 
any either. Indeed, it turns out that the maximum in Eq. (|60|) is assumed for (62 = 0, 63 = 0). 
Therefore we will only work with the nontrivial variational parameters 69 and bi. Because 
the variational eigenfunction is linear in the parameters, there is an equivalent method that 
we will use. We construct an orthonormal basis out of 1, \v\'^ and \v\^i and determine the 
matrix elements of the operators L and W^^^. The largest eigenvalue of the truncated 



matrix (in which all other matrix elements are set to zero) 
bound to the real eigenvalue Aq. 
The basis vectors are 



is an approximation and lower 



|1) = 1, |2) 



|3) 



+ 1, 



which are orthonormal with respect to the inner product defined in Eq. 
straightforward calculation gives the truncated matrices on this basis. 



). A tedious but 
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FIG. 3. Several approximations for the function w. A 1 x 1 truncation of the matrices yields the 
dotted curve. The solid curve results from the 3x3 matrices. As the difference is so small, only 
the blow-up in the inset can show the 2x2 results separately as the curve of points. 

For many values of Wq'j the eigenvalue problem 

was solved numerically using "Maple V", and the largest eigenvalue, Ao(wo7) was taken to 
determine w, according to Eq. (|62D . The minimum wq is assumed at about tfoTo = 3.5, 
therefore, in Fig. ^, the values of w were plotted for 2 < < 4. 
The minimum of w can be determined numerically. It is given by 



Wq ^ 4.735 . . . , 
woio ~ 3.497..., 

and the normalized eigenfunction is approximately 

Aq{v) = 0.63943 + 0.16289|t;|2 + 0.004351|w|l 



(64) 



(65) 



In Fig. ^, also the result is plotted when the matrices are truncated further. Keeping only |1) 
means no variational parameters and no velocity dependence. Hence we recover the result 
Wo ~ 4.31107. . .of previous work [^,|l^], in which the velocity dependence of the collision 
frequency was neglected. The enhancement of taking two basis vectors, |1) and |2), is rather 
large, but taking one more gives only a small difference, of about 0.1%, as the inset of Fig. |^ 
shows. Therefore we assume that these values have converged up to at most 0.1%. 

The results can be qualitatively understood as follows. In the first place, particles with 
a higher velocity have a higher collision frequency, and so their clock values increase faster 
than average. Thus, one expects higher velocities to be more prominent in the head than in 
the rest of the distribution. Indeed, this is seen in Eq. (§^), where according to formula (|53[), 
the positive value of the coefficient 0.16289 signifies a shift to higher velocities in Phead as 
compared to ipo. Secondly, collisions of other particles with the head just synchronize this 



24 



particle to the one in the head. So the colhsion frequency of the particle in the bulk is 
irrelevant and unable to compensate the increase in collision frequency in the head. This 
increase is what really causes the increase of w over the velocity independent estimate of 
4.31107. . . in Refs. [jl5|,|l6|, as can be seen from the relative increase in this frequency, 



^head 
^d 



20F 



1 



1. 15^ 
2^0 + 



1 + -&n + — &i = 1.0977..., 
which matches closely the relative increase in clock speed, wo/4. 31107 = 1.0983 . . .. 



B. Comparison with Simulations 

The results for the largest Lyapunov exponent of a system of hard balls will now be 
checked against simulations. We want to check the following points: the validity of the clock 
model, the relation between the Lyapunov exponent and the clock speed, and the velocity 
distribution in the head. 

The simulation method we shall use is a variant of the Direct Simulation Monte Carlo 
Method (DSMC), that was also used by Dellago and Posch [^] for the calculation of Lya- 



punov exponents in a "spatially homogeneous system". These authors also checked the 
equivalence for low densities of DSMC and molecular dynamics simulations. The low densi- 
ties that we are interested in, are not accessible to the molecular dynamics simulations. The 
DSMC method simulates the Boltzmann equation, and consists of the following. In each 
time step, a pair of particles is picked at random to collide. The probability of accepting a 
pair is proportional to the relative velocity so that the collision frequency of two particles 
with given velocity is also linear with the relative velocity. In the original method due to 



Bird the system is divided into cells, and only particles within one cell can collide, but 
as we are simulating a homogeneous system, the whole system is in one cell. Once a pair is 
picked, a collision normal a is drawn from a distribution proportional to \v2\ ■ (3"|, like in the 
Boltzmann equation. The velocities of the pair are transformed according to Eq. (|^) and 
their deviations are subjected to the transformations in Eq. (0), to a account for free flight, 
and, subsequently Eq. (|^). Even though in the DSMC method the positions themselves are 
discarded, we still keep track of the deviations in position by integrating the deviations in 
velocity by means of Eq. (^). 

In a parallel simulation, for the same collision sequence, the particles are given separate 



clock values, which are updated in a collision according to Eq. (|5lD . This enables us to check 
how well the clock values from the dynamics of Eq. (^) represent the actual dynamics of 
the velocity deviations. 

Initially, the velocities of the particles are picked from a Maxwellian distribution, scaled 
by Vo, i.e., with /c^T/m set to one. The initial velocity and position deviations are unit 
vectors drawn from an isotropic distribution, 5fo is set to one (this can be done because the 
deviations follow a linear equation), and correspondingly all clock values are zero. 

The first check is to see whether the clock values are accurate in describing the behavior 
of the velocity deviations. As mentioned above, in the simulation the particles have their 
real deviations in position and velocity as well as a clock value, evolving independently of 
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the deviations, according to Eq. (^). If that equation were exact, for each particle the 
clock value ki would equal In \6vi/6vo\/\ lnn| for all time, if it did so initially. But Eq. (|50D 
shows that there are density corrections of order 1/ In n. If the difference between and 
In |(5i7j/(5uo|/| lnn| became very big too often, the maximum of the clock values might not 
correspond to the maximum of the velocity deviations, and the dynamics from Eq. ( |5TD 
would not even be approximately correct. If things go as we expected, the clock values will 
give the right w in the limit of zero density, i.e., wq, with deviations that scale as 1/lnfi. 

Two simulations for N = 128 particles, and d = 2 dimensions, were performed for, in 
total, 3, 887, 196 collisions, one simulation at a density n2 = 10~^ and one at a density 
77,2 = 10~^^. In Fig. ^ results are plotted for both simulations. For each particle a symbol is 
plotted, a plus for the simulation at density n2 = 10~^^ and a diamond for the simulation 
at density n2 = 10~^. The horizontal position of each symbol is given by the final clock 
value ki of particle i (of which 226885 is subtracted), the vertical position by the final value 
of In l^-Uj/^fol/l lnn2|. The scale of the velocity deviations is a bit different for the two 
simulations, therefore different scales were used on the vertical indicated in the 

figure. 

From the plot, we can see the following. There is a linear relation between 
In \6vi/6vo\/\ In 77,2! and fcj, around which there are fluctuations. The fluctuations get smaller, 
and the slope gets closer to one, as the density decreases. Small fluctuations indicate that 
the maximum in Eq. (|5TD is correct, i.e., the particle with the largest clock value has the 
largest velocity deviation. The value of In \ 6vi/6vQ\/\ lnn2| is not exactly equal to the clock 
value ki. Eq. (|50D shows that the correction in the increase of the clock value in a colli- 
sion is of the order of l/lnn2. Because such a correction is present in each collision, the 
correction to the relation In \ 6vi/SvQ\/\ lnn2| = ki also grows. The relative correction is, for 
density ^2 = 10~^ of the order of 24%, while for the density n2 = 10^^^ it is about 8%. 
The difference by a factor of three shows that the correction is of the order of 1/| lnn2|, as 
In 10~^^/ In 10~^ = 3 also. We conclude that the dynamics of the clock value in Eq. ( |^ 
indeed describes the dynamics of the In \ Svi\ accurately for low densities. 
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FIG. 4. Check on the clock model, for N = 128 and d = 2. For each particle, |(5^i|/| lnn2| is 
plotted against /cj. The diamonds are from a simulation with n2 = 10~^, after 3, 887, 196 collisions. 
The plusses are from a simulation with 77-2 = 10~^^. Also plotted is a line with slope 1. 



The second check is on the asymptotic form (^) of the Lyapunov exponent, and its 
relation to the leading clock speed. To check the leading behavior of the largest Lyapunov 
exponent as a function of density, we need to perform a number of simulations for different 
densities, all very low, because the coefficient in the next term in Eq. (0), although it could 
be calculated in principle, is not yet known and thus we cannot assume that we can neglect 
it. The leading clock speed wq, which can be determined from the parallel calculation of 
the clock values, Eq. (|5T|) , does not depend on density and can be determined from just one 
simulation at an arbitrary density. 

In this case, simulations are done with N = 64 particles, again in d = 2 dimensions. The 
Lyapunov exponent divided by the collision frequency is plotted as a function of the density 
77-2 in Fig. ^. To justify the form of the density expansion given in (0), we also plotted the 
function —wolnn2 + Wi, where wq is the clock speed found from the parallel simulation of 
the clock values, Wq = 3.4479 ± 0.0016 and Wi is obtained from fitting (using only points 
from densities 77,2 smaller or equal to 10~^). This function describes the simulation points 
very well, i.e., the slope is indeed given by the clock speed. The deviations at higher density 
may be accounted for by higher powers of 1/ lnn2. 

Surprisingly, the value of wq is far from the predicted value in Eq. (p^. This is due to 
finite size effects, which have been well studied by now p!6| , p9| , |30| . Corrections due to the 
finiteness of N are, for large N, of the order of ln~^ A^, but for iV = 64 we are not yet in 
the asymptotic regime for which this scaling holds. The investigation for larger numbers of 
particles and the calculation of the finite size corrections will be done in future work, but 
we remark that for a simplified clock model, in which the collision frequency is taken to be 
velocity independent, this was already carried out in 
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FIG. 5. Density dependence of the largest Lyapunov exponent A,: 



0.01 



for iV = 64, d = 2. The 



plusses are DSMC results. The linear curves is a function of the form (52), where the first two 
terms are included. For wq, the clock speed from the simulation is taken, wi is a fit parameter. 
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Thirdly, we will check the velocity distribution Phcad(^^) in the head of the clock distri- 
bution. This distribution is determined as follows. After the simulation has run for some 
time, we determine the average of the clock values and shift all the clock values by that 
amount, so we can access the stationary F{k, v) rather than the shifting C{k, v, t). We then 
pick some lowest clock value fco; all particles with a lower clock value are discarded. Of the 
remaining particles, which are those in the head of the distribution, the velocity distribu- 
tion is constructed. We adjust ko such that k > ko can be identified with the tail of the 
distribution. To get good statistics (which is difficult because there are not many particles 
in the head), we measure this distribution at several times, and average the result. 

As the velocity distribution Phead{v) only involves clock values, there is no density de- 
pendence and a simulation at one density is enough. We ran a simulation with N = 128 
particles, and d = 2. The distribution of the clock values was obtained, and a value of fco = 7 
seemed a reasonable start of the tail of the distribution. We have checked that the results 
do not change much when instead we take k^ = Q or ko = 8. 

Combining Eqs. (|65|) and (|63| ) gives the explicit prediction for the velocity distribution 
in the head: 

Phcad(^) = (0.63943 + 0.16289t;2 + 0.004351^;^) v e-^^^^. (66) 

This prediction has been plotted in Fig. |] together with the distribution found from the 
simulations, and the velocity distribution in the whole gas. Even though the statistics isn't 
perfect, there is a very good agreement between the two. It should be noted that leaving 
out the term in Eq. (p6| ) gives only a small difference. It is also evident that the velocity 
distribution in the head of the clock distribution has higher mean velocity than the velocity 
distribution characterizing the full gas, which is also plotted in figure Eq. (||). 




FIG. 6. Comparison of the velocity distribution in the leading edge. The histogram results from 
the simulations for N = 128, in d = 2 dimensions. The solid curve is a plot of the prediction 
The dotted curve is the two-dimensional Maxwell distribution, t; exp(— z;^/2). 
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V. THE DILUTE, RANDOM HARD-BALL LORENTZ GAS 



A very useful and simple example of a dynamical system which is a dispersing billiard, 
rather than a semi-dispersing billiard is provided by the Lorentz gas model. In this model 
one places fixed, hard balls of radius a in a d-dimensional space, and then considers the 
motion of a point particle of mass m and with initial velocity v, in the free space between 
the scatterers [Q. The particle moves freely between collisions and makes specular, energy 
conserving collisions with the scatterers . We consider here only the case where the scatterers 
are not allowed to overlap each other. One can imagine that the scatterers are placed with 
their centers on the sites of a regular lattice, or that the scatterers are placed randomly 
in space. Here we consider the case where scatterers are placed randomly in space in a 
volume V in such a way that the average distance between scatterers is much greater than 
their radii, a. That is, the scatterers form a quenched, dilute gas with na'^ -C 1. 

In this section we show that the simplifications in the dynamics of this model produce 
simplifications in the calculations of the Lyapunov exponents and KS entropies of the moving 
particle, due primarily to the dispersing nature of the collisions of the moving particle with 
the scatterers. A Lorentz gas in d dimensions may have at most d — 1 positive Lyapunov 
exponents. That is, the phase space dimension for the moving particle is 2d, the requirement 
of constant energy removes one degree of freedom, and the direction in phase space has a 
zero Lyapunov exponent associated with it, since two points on the same trajectory will not 
separate or approach each other in the course of time. Thus there can be at most 2d — 2 
nonzero Lyapunov exponents, of which only d — 1 can be positive, since the exponents come 
in plus-minus pairs for this Hamiltonian system. 

To analyze the Lyapunov spectrum we consider two infinitesimally close trajectories in 
phase space and follow the spatial and velocity deviation vectors 6r, 5v separating the two 
trajectories, in time. By arguments almost identical to (but simpler than) those presented 
in Section II for the case where all particles move, we have the following equations of motion 
for the position, r, velocity, v spatial deviation vector, 5r, and velocity deviation vector, 5v, 
of the moving particle, between collisions 

r = V, 
6r = 6v, 

6ij = 0. (67) 
At a collision with a scatterer, these dynamical quantities change according to 

f, 

V — 2{v ■ o")(3", 
■ 6f, 

M^-Sv- 2Q<j ■ Sr. (68) 

Again, the primed variables denote values immediately after a collision, while the unprimed 
ones denote values immediately before the collision. The distance from the center of the 
scatterer to the moving particle, at collision, is aa, and the tensor is given by 



r 
v' 

6v 
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Q_ _ [{b-v)l + av\ ■ [{a-v)l-va] 
a{cr ■ v) 

It should be noted that if the velocity deviation vector and the velocity vector are orthogonal 
before collision, then they will be orthogonal after collision as well. That is, if w ■ 5w = 0, 
then if -Stf = 0, also. Thus we can take 6v to be perpendicular to v for all time, and without 
loss of generality, we can take 6r to be perpendicular to v as well. Of course, the condition 
that V ■ 6v = is simply the statement that the two trajectories are on the same constant 
energy surface. We will denote the spatial and velocity deviation vectors for the Lorentz 
gas with a subscript _L to indicate that they are defined in a plane perpendicular^ to v. It 
follows that we may replace the d x d ROC matrix defined by 

6f=p-6v (70) 

by a ((i — 1) X (d — 1) matrix p_|_ defined by 

6f± = p^-5v±. (71) 

For the two dimensional case pj^ is a simple scalar which we denote by p. One easily finds 
that between collisions p grows with time as 

p{t)=p{0)+t, (72) 



The change in p at collision satisfies the "mirror" equation |T8 



1 1 2t; 
p' p a cos (p 

Here the primes denote values immediately after a collision, and v is the magnitude of the 
velocity of the particle. If the radius of curvature p is positive, initially, it will always remain 
positive, and it also follows from Eq. ( [73| ) that the value of vp after collision is less than half 
of the radius of the scatterers. Consequently the radius of curvature typically grows to be 
of the order of a mean free time, and it becomes much smaller immediately after a collision 
with a scatterer. For three dimensional systems a similar situation results. Now p_|_ is a 
2x2 matrix which satisfies the free motion equation 

pAt) = P±{0)+tl, (74) 

and changes at collision according to 




va + av — — — — — [v ■ (Til 



■ M^. (75) 



[V ■ a) 

Here the inverse radius of curvature matrices [p^]^^ ^^'^ defined in planes per- 

pendicular to v' and to v, respectively. In the hybrid notation of Eq. ([75|) , in which both 



^Of course, the components of 5v, 6 fin the direction of v are not related to the non-zero Lyapunov 
exponents or the KS entropy, since these components do not grow exponentially. 
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dx d matrices and d — 1 x d — 1 matrices figure, the inverse matrices in the directions along 
v' and V, respectively, may be defined by [p'_^_]~^v' = and [p^]"^^ = 0. The final matrix 
in the right-hand side of Eq. ( [75| ) can then be restriced to the plane perpendicular to v' 
straightforwardly. 

It is worth pointing out some important differences between the ROC matrices defined 
here for the Lorentz gas and those defined earlier for the regular gas of moving particles. 
Here the ROC matrices are defined in a subspace orthogonal to the velocity of the moving 
particle. Further the change in the matrix elements at collision is from a typically large 
value on the order of a mean free time to an always small value, on the order of the time 
it takes to move a distance equal to half the radius of a scatterer. This latter property is a 
property of dispersing billiards. For the regular g cLS CcLSG, only a few of the elements of the 
ROC matrices become small after a collision, which means that one cannot find an accurate 
approximation to the ROC matrices by considering only one collision. This latter property is 
associated with semi- dispersing billiards where a refiection from a scatterer does not change 
the diagonal components of the ROC matrix that correspond to the fiat directions of the 
scatterer, at all. 



A. Informal Calculation of the KS Entropy and Lyapunov Exponents for the Dilute, 

Random Lorentz Gas 

Here we show that simple kinetic theory methods allow us to compute the Lyapunov 
exponents and KS entropies of Lorentz gases in two and three dimensions . To do that 
we use methods similar to those in Section HI. That is, we consider the equations for the 
deviation vectors, Eqs. (^7\ - ^U]) above. The velocity deviation vector changes only upon 



collision with a scatterer. We will base our calculation on the exponential growth rate of 
the magnitude of the velocity deviation vector, and for three dimensional systems, on the 
exponential growth rate of the volume element in velocity space. 

We begin by writing the spatial deviation vector just before collision as 

6f=t6v + 6f{0), (76) 

where 6r{0) is the spatial deviation vector just after the previous collision with a scatterer. 
This equation is essentially the same as Eq. (|1^), but now it is a good approximation to 
neglect the spatial deviation vector vector 6f{0) since 6f{0) is of relative order a/vt compared 
to the term t6v, in all directions of 6f, perpendicular to v. Thus we neglect this term and 
insert Eq. ( [76D into the last equality of Eq. (|68D to obtainQ 

dv' = Ma- 6v - 2tQa ■6v = Si- 6v, (77) 

where we have defined a matrix a that gives the change in the velocity deviation vector at 
collision. Then we can express the velocity deviation vector at some time t in terms of its 
initial value as 



^In principle, the term M5- • 6v in Eq. ( [77| ) can be neglected also, but only for directions perpendic- 
ular to the velocity. If one is careful to consider only deviations 6v in the subspace perpendicular 
to the velocity of the particle, it is possible to carry out the calculation with this term neglected. 
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Sv(t) = a^f ■ a^r-i ■ ■ ■ ai • Sv{0), (78) 

where we have labeled the successive collisions by the subscripts 1,2,..., A/". We can de- 
termine the largest Lyapunov exponent by examining the growth of the magnitude of the 
velocity deviation vector with time, and the KS entropy as the growth of the volume element 
with time. Therefore, with the approximations mentioned above, 

Am«.r = lini - In 



= hm — -7 V In ,' J; ', , (79) 

where 5vf is the velocity deviation vector immediately after the collision labeled by the 
subscript i. Similarly, the sum of the positive Lyapunov exponents is given by 

U 1 ^ 



Ai>0 



t 



To evaluate the sums appearing in Eqs. (|79| , ^), we note that to leading order in the 



density none of the collisions are correlated with any previous collision, that is, the leading 
contribution to the Lyapunov exponents comes from collision sequences where the moving 
particle does not encounter the same scatterer more than once in the sequence 0. Therefore 
we can treat each term in the sums in Eqs. (^, ^D|) as being independent of the other 



terms in the sum. We have expressed Xmax and the sum of the positive Lyapunov exponents 
as arithmetic averages, but for long times and with independently distributed terms in 
the average, we can replace the arithmetic averages by ensemble averages over a suitable 
equilibrium ensemble. That is 



and 



\5v 



(81) 



^ Ai = z/(ln|deta|) (82) 

Ai>0 

where 5v~ and 5v^ are the velocity deviation vectors before and after collision, respectively, 
V is the (low density)value of the collision frequency, M /t as t becomes large, and the angular 
brackets denote an equilibrium average. 

We now consider a typical collision of the moving particle with one of the scatterers. The 
free time between one collision and the next is sampled from the normalized equilibrium 
distribution of free times P|, -P(t) given at low densities by 



®In two dimensions the particle will hit the same scatterer an infinite number of times. However 
the effects of such processes are of higher density, and can be neglected here since the number of 
collisions between successive collisions with the same scatterer become typically very large as the 
density of scatterers approaches zero. 
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P(t) 



(83) 



The construction of the matrix a requires some geometry and depends on the number of 
dimensions of the system. In any case we take the velocity vector before colhsion, v to be 
directed along the 2;-axis, and take a ■ v = — f cos0, where — 7r/2 < < 7r/2. The velocity 
deviation before collision 6v~ is perpendicular to the 2;-axis. Then it is a simple matter 
to compute and | deta|. For two-dimensional systems 6v and the matrix a are 

given in this representation hyf\ 



6v' 




1 + A) cos 20 sin 20 
[1 + A) sin 20 -cos 20 



(84) 



where we introduced A = (2fr)/(acos0). To leading order in vr/a we find that 



A; 



det a| = A. 



(85) 



For three dimensional systems the unit vector a can be represented as a = — cos 5 + 
sin0cosa x + sin sin a y. Now the ranges of the angles and a are < < 7r/2 and 
< a < 2tt. There is an additional angle ip in the x, y plane such that the velocity deviation 
before collision 6v~ = \ 6v~ \ [x cos ip + y sin ip] . It is somewhat more convenient to use a 
symmetric matrix, a = (1 — 2aa) ■ a, given by 



/ 1 + A(cos^ + sin^ cos^ a) 

A sin^ cos a sin a 
\ 



A sin^ cos a sin a 
1 + A(cos^ + sin^ sin^ a) 
1 



(86) 



One easily finds 



\Sv- 



2tv 



\Sv- 



cos^(a — ip) 



cos^ 



+ sin^(Q; — -0) cos^ 



1/2 



and 



det a| = I det a| 



(87) 



(88) 



to leading order in vr/a. 

To complete the calculation we must evaluate the averages appearing in Eqs. (|7^, pOl ). 
That is we average over the distribution of free times and over the rate at which scattering 
events are taking place with the various scattering angles. Additionally in 3 dimensions an 
average over a stationary distribution of angles ip has to be performed in general. Due to 
the isotropy of the scattering geometry ip can here be absorbed in a redefinition a' = a — ip 



^In contrast to the ROC matrices p, a is a d x d matrix. If one chooses one of the basis vectors of 
a perpendicular to v, the remaining ones are the basis of the corresponding d — 1 x d — 1 matrix, 
from which one can also obtain Eq. (^5|), and, in the three dimensional case, Eqs. ( |87|) and (|88|). 
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of the azimuthal angle a. This will not be true any more if the isotropy of velocity space is 
broken (e.g. by an external field). The appropriate average of a quantity F takes the simple 
form 

(F) = - r dr [ da cos (pP{T)F, (89) 
J Jo J 

where -P(t) is the free time distribution given by Eq. ( p3D and J is a normalization factor 
obtained by setting F = 1 in the numerator. The integration over the unit vector a, i.e., 
over the appropriate solid angle, ranges over — 7r/2 < < 7r/2 in two dimensions and 
over < < 7r/2 and < a < 27r in three dimensions. After carrying out the required 
integrations we find that 

= Xmax = 2nav[- \n{2na^) + 1 - C] + ■ ■ ■ , (90) 

for two dimensions. Here C is Euler's constant, and the terms not given explicitly in Eq. 
(pop are higher order in the density. Similarly, for the three dimensional Lorentz gas we 
obtain 

A+a. = na'v7r[- ln(n/2) + In 2 - 1 - C] + ■ ■ ■ , (91) 



ALx + A+,„ = 2na\7f[- \n{n/2) -€] + ■■■, (92) 
from which it follows that 

A+,„ = na^nl- ln(n/2) - In 2 + ^ - C] + ■ • • , (93) 
where n = na^ir. We have therefore determined the Lyapunov spectrum for the equilib- 



rium Lorentz gas at low densities in both two and three dimensions |T^,[T^. There is good 
agreement with simulations, as shown in figure |^. We note that the two positive Lyapunov 
exponents for three dimensions differ slightly, and that we were able to get individual values 
because we could calculate the largest exponent and the sum of the two exponents. We 
could not determine all of the Lyapunov exponents for a d > 3 dimensional Lorentz gas 
this way. Moreover, for a spatially inhomogeneous system, such as those considered in the 
application of escape-rate methods, the simple kinetic arguments used here are not sufficient 
and Boltzmann-type methods are essential for the determination of the Lyapunov exponents 
and KS entropies. 

In Fig. (|^) we illustrate the results obtained above for the Lyapunov exponents of the di- 
lute Lorentz gas in both two and three dimensions, as functions of the reduced density of the 
scatterers, and compare them with the numerical simulations of Dellago and Posch pT| , |52 
As one can see the agreement is excellent. 
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FIG. 7. A plot of the Lyapunov exponents, in units of v/a, for the moving particle in a random, 
dilute Lorentz gas in two dimensions (top) and three dimensions (bottom), as functions of the 
density n, in units of a"'^. The solid lines are the results given by kinetic theory, Eq. (pO|), 
respectively, Eqs. (pl|) and (|93|), and the data points are the numerical results of Dellago and 
Posch. 

B. Formal Kinetic Theory for the Low Density Lorentz Gas 

The formal theory for the KS entropy of the regular gas is easily applied to the Lorentz 
gas, which is, of course, considerably simpler. Thus by following the arguments leading to 
Eq. ([T3|) for the sum of the positive Lyapunov exponents for the regular gas, we find that 
the KS entropy for the equilibrium Lorentz gas is given by 

^ Ai = a'^^^ / dx dpdRdaQ{—v ■ a)\v ■ a|(5(r — R — aa) x 

A,>0 
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In I det [M^ + 2Q^ ■ p] |JS(2;, R, p). 



(94) 



Here R denotes the location of the scatterer with which the moving particle is colliding, and 
T2 is the pair distribution function for the moving particle to have coordinate, r, velocity, 
w, ROC matrix, p, while the center of the scatterer is located at R. At low densities we 
may assume that the moving particle and the scatterer are uncorrelated, so that the density 
expansion for JF2, immediately before collision would have the form 



J^2{x,R,p) = nTi{x,p) 



(95) 



where n is the number density of the scatterers and J^i (x, p) is the equilibrium single particle 
distribution function for the moving particle. 

We may easily construct an extended Lorentz-Boltzmann equation (ELBE) for along 
the lines used previously for the extended Boltzmann equation in Section III. The ELBE is 
given by 



d_ 

dt 



J^i{x, p) 



dRda&{—v ■ cr)\{v ■ a) x 



6{f- R-aa) / dp'5{p - p{p'))V'^ -5{f-R + aa) 



Ti{x,p). 



(96) 



Here the operator is a substitution operator that replaces velocities, v and ROC matrices, 
p, by their restitution values, i.e., the values needed before collision to produce the values 
and p after collision with a scatterer with collision vector a. Also, the free particle streaming 
operator on the left-hand side of Eq. (|96|) is given by 



d 



— + y 

df ^1 dp 



(97) 



since between collisions, r varies as r(0) + tv and p varies as p(0) + tl. 

Returning for the moment to Eq. ( p^ ) for the KS entropy, we see by evaluating the 
determinant in the integrand on the right side, that not all components of p contribute. In 
fact, only the components of p in the plane perpendicular to v contribute to hj^s- Here we 
will work out the details of the calculation of hxs for the two dimensional case, leaving the 
details of the three dimensional case to the literature WM. For d = 2 we can easily evaluate 



the determinant in Eq. (|94D, and find that it is 

I det [M^ + 2Qa-p]\ = 



2vp 



acosi 



(98) 



Here the scalar p is the component of the p matrix given by p = -Oj, ■ p ■ -0^, where -0^ is 
a unit vector orthogonal to v. Moreover, for low densities we can use the approximation 
nJ^i{x, p) for J^2- Then the expression for h^s becomes, at low densities 



hRS = an J dxdpdaQ{—v ■ a)\v ■ a\ In 



2vp 
a cose 



Ti{x,p), 



(99) 



36 



where x = f,v, and we have to determine Fi{x,p) as the solution of the ELBE where we 
integrate over all components of p, except the one diagonal component p. The ELBE then 
becomes, in the spatially homogeneous, equilibrium case 



d_ 

dp 



J^i{f,v,p) 



nav 



Tr/2 

I 

-7r/2 



COS( 



acos( 



"^'T 2v + - 



\j^i{f,v',p') -J='i{f,v,p) 



vp' 



(100) 



The argument of the 6 function is simply obtained by using the mirror formula given by Eq. 
([f3|), but now the unprimed variable is the value a/ier collision, and the primed variable is the 
value of p before collision. A further, and useful simplification results from the observation 
that p' is typically of the order of the mean free time between collisions, which, for low 
density, is much larger than a/v. Therefore the delta function can be replaced by 



6 p 



a cos I 



2v + 



6 p- 



acos( 



vp' 



2v 



(101) 



In a spatially homogeneous and isotropic equilibrium state, J-'i has to be a function of the 
norm of the velocity vector \v\ and the radius of curvature p only. Using that the magnitude 
of the velocity, \v\, always stays the same, we know that there is a solution of the form 



^i(,r,v,p) = i^{v)i){p), 



(102) 



with ^{v) = {27iV)~^5{\v\ — Vq) is the normalized equilibrium spatial and velocity distribu- 
tion function for the moving particle, Vq is its constant speed, and V is the volume of the 
system. Now all we have to do is to determine ip{p). An inspection of Eq. ( p.OU| ), with the 
approximation, Eq. ( |101|) shows that ip{p) satisfies the equation 



dp 



+ 2navip{p) = 0, 



for p > a/{2v), with solution 

V'(P) 



(l/to)e~*/*° for p > a/{2i 



(103) 



(104) 



Here to is the mean free time given by to = {2nav) ^. In the case that p < a/{2v), one can 
easily solve the full equation with the delta function to find 



^(p) = (1/to) 1 



2vp 



1/2 ~ 



for p < a/ (2t>). 



(105) 



Combining these results with Eq. ( |102| ) and inserting them in Eq. ( p9|) we recover the 
result, Eq. ( ^0]) for the KS entropy, equivalently the positive Lyapunov exponent, for the 
low density, equilibrium, random Lorentz gas. A similar, but somewhat more elaborate 
calculation can be carried out, for the three dimensional case, to obtain exactly the same 
result as obtained in the informal theory for the KS entropy. To obtain the largest Lyapunov 
exponent by more formal methods, one has to resort to methods for determining the largest 



eigenvalue for products of random matrices. This is well described in the literature and 
we will not pursue this issue further here. 
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VI. CONCLUSIONS AND OPEN PROBLEMS 



In this article we have given a survey of the apphcations of the kinetic theory of dilute, 
hard-ball gases to the calculation of quantities that characterize the chaotic behavior of such 
systems. Results have been obtained for the KS entropy per particle and for the largest 
Lyapunov exponent of dilute hard-ball systems, as well as for the Lyapunov spectrum for 
the moving particle in the dilute, random Lorentz gas, with nonoverlapping, fixed, hard- 
ball scatterers. All of these results are in good to excellent agreement with the results of 
computer simulations. In the study of the largest Lyapunov exponent, we have developed 
a very interesting clock model which seems to explain many features of the behavior of 
this exponent. Moreover, the method for treating the clock model reveals a deep and, 
perhaps unexpected, connection between the theory of Lyapunov exponents and the theory 
of hydrodynamic fronts. 

We emphasize that the results given here apply to a dilute gas in equilibrium, but 
nonequilibrium situations have been treated by these methods as well. For example, it is 
possible to calculate the Lyapunov spectrum for a dilute, random Lorentz gas in a nonequi- 
librium steady state produced by a thermostatted electric field, at least for small fields, and 



to obtain results that are in excellent agreement with computer simulations [^,0,0 . Cal- 
culations are currently underway for the largest Lyapunov exponent for a hard-ball gas, when 
the gas is subjected to a thermostatted, external force that maintains a steady shear flow 
in the gas. Furthermore, one can use kinetic theory to calculate the Lyapunov spectrum for 



trajectories on the fractal repeller for a Lorentz gas with open, absorbing boundaries ||13| , |35 
Such results are useful for understanding escape-rate methods for relating chaotic quanti- 
ties for trajectories on a fractal repeller of an open system to the transport properties of 
the system, as described by Gaspard and Nicolis P, p!0| , p6| . These results will eventually be 
extended to hard-ball gases as well. 

Of course, many problems remain to be solved. Here we mention some of the most 
immediate problems: 

1. All of the results described here apply to hard-core systems. That is the particles 
interact with a potential energy that is either zero, beyond a given separation, or 
infinite, below that separation. It is worth studying the properties of dilute systems 
with smoother potential energies for a number of reasons: (a) The results of Rom- 
Kedar and Turaev suggest that the dynamics of particles with short ranged but 
smooth potentials may exhibit regions of nonhyperbolic behavior. It is important to 
know more about these regions and to assess their effect on the overall chaotic behavior 
of gases interacting with such potentials, (b) We know very little about the chaotic 
properties of dilute gases interacting with long range forces, such as Maxwell molecules 
and Coulomb gases. 

2. An open problem, even for dilute gases, is to obtain the complete spectrum of Lya- 
punov exponents for a hard-ball gas. This is a very challenging problem in mathe- 
matical physics, and no easy approach is in sight. There is a very tantalizing set of 



numerical results by Posch and coworkers ||3^ showing that the smallest nonzero Lya- 



punov exponents have a hydrodynamic structure, in that the exponents themselves 
seem to scale as the inverse of the linear size of the system, and that the spatial 



38 



deviation vectors seem to form collective modes of both transverse and longitudinal 



types. Eckmann and Gat |^ have proposed an explanation of this hydrodynamic-like 
behavior of the lowest Lyapunov modes using techniques from the theory of random 
matrices. It would be useful to understand and to extend their results using kinetic 
theory methods. 

3. The extension of these results to gases at higher densities remains an open, challenging 
problem. The kinetic theory of dense gases has exposed a number of effects caused 
by long range dynamical correlations between the particles. Such effects include non- 
analytic terms in the density expansion of transport coefficients and long time tail 
phenomena in time correlation functions, among others 0]. It is of some interest to 
see the effect of these dynamical correlations on the chaotic properties of the gas, as 
well. Furthermore, a high density hard-ball system may form a glass, and there would 
be useful information obtained about the glassy state if one could study the chaotic 
behavior of the gas through the glass transition. 

4. The extension of the clock model to higher density systems can also be expected to 
reveal new and interesting phenomena connected to the effects of density and other 
fluctuations in the gas upon the clock speed and related quantities. At present there 
are some weak numerical indications that the largest Lyapunov exponent might 



diverge in the limit of large numbers of particles as InA^, where N is the number of 
particles in the system. It might be possible to confirm or to rule out this possibility 
by extending the clock model so as to include the effects of fluctuations in the fluid. 

5. Lyapunov exponents and KS entropies are not the only properties that characterize 
chaotic systems. There are many more quantities, such as topological pressures, fractal 
dimensions, etc., that remain to be explored by the methods outlined here. 

6. A subject of considerable interest and activity is the physics of gases that make inelastic 
collisions, i.e. the physics of granular materials |^0[. One would like to know what the 
chaotic properties of such systems might be, or, more generally, how to deflne such 
properties in nonstationary systems. 

In conclusion, we have described in this article only the flrst steps, taken over the last 
few years, to develop useful methods for the calculations of Lyapunov exponents and KS 
entropies for the systems of particles that can be treated by kinetic theory. We are partic- 
ularly delighted that kinetic theory has something to contribute to the fleld of the chaotic 
behavior of large systems of particles, and that the ideas of Maxwell and Boltzmann are still 
having new and fruitful applications. 
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